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Abstract 


In polymer electrolyte fuel cell (PEFC), it is important to understand the behavior of liquid water in gas diffusion layer (GDL) which is one 
of the constructional elements so as to improve the output performance and the durability. As this behavior of liquid water is attributed to not 
only the hydrophilicity but also inhomogeneous structure, it is needed to examine in consideration of an actual GDL structure. In this study, 
as the basic examination of two-phase flow analysis in an actual GDL, a simulated GDL was made by numerical analysis considering the fiber 
placement. Furthermore, the prediction methods for pore size distribution, permeability and tortuosity of this simulated GDL were developed with 
the numerical analysis. These parameters of flow and mass transfer were compared with other studies, and the validity of this simulated GDL 
was confirmed. In addition, effective diffusion coefficient was calculated from tortuosity in simulated GDL, and PEFC output performance was 
evaluated by a simple model. Moreover, the optimal GDL was examined in consideration of the effect of porosity and fiber diameter at the fiber 


level. 
© 2007 Elsevier B.V. All rights reserved. 


Keywords: PEFC; Numerical analysis; Gas diffusion layer; Pore distribution; Permeability; Tortuosity 


1. Introduction 


Recently, global energy consumption increases, and global 
environmental pollution and destruction of ecosystem turn 
worse by mass consumption of fossil fuels such as petroleum. 
The exhaustion of these energy resources becomes a serious 
problem. In order to solve these problems, the practical applica- 
tion of fuel cells are expected because it emits less environmental 
pollutant and converts more efficiently from chemical energy to 
electrical energy than other energy resources. There are various 
kinds of fuel cells, and polymer electrolyte fuel cell (PEFC) is 
especially expected as driving power of vehicles and stationary 
power supply, because it can work at low temperature and has 
high power density. The performance of PEFC has improved 
rapidly by developing the new component materials and opti- 
mizing the system, and some PEFC systems have already been 
commercialized. However, in order to spread PEFC for various 
uses, it is necessary to improve the durability and the cell out- 
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put, and to reduce the cost. In addition, it is needed to investigate 
the phenomena which influence the output performance and the 
durability of PEFC. The power generated by PEFC is affected by 
the structure, the material and the operating conditions through 
the process of generation. The phenomena of mass transfer, heat 
transfer, catalysis, electrochemical reactions and fluid dynamics 
are shown only in an internal cell, and it is greatly important 
to understand the correlation among such complex phenomena 
in detail to improve and optimize the PEFC component and the 
system. 

The phenomena and the mechanism in gas diffusion layer 
(GDL) which is one of the PEFC components have been studied. 
Theroles of GDL are the followings: the electro transfer between 
separator and catalyst layer; the mass transfer of hydrogen, oxy- 
gen and water between gas channel and catalyst layer; enlarging 
reaction area of electrode. Accordingly, GDL should be porous 
media and have high electrical conductivity. In general, GDL is 
made from carbon materials; carbon paper (non-woven fabric) 
and carbon cloth. They are distinguished by the difference of the 
laminated structure of a carbon fiber which is several microme- 
ters in diameter. The improvement of mass transfer rate in GDL 
is important to generate high power density, and problems of 
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Nomenclature 


Cm 
channel 


o 
S 


WT PUR li 


velocity of virtual particle moving 

oxygen concentration at a channel (mol m^?) 
oxygen concentration at an interface of catalyst 
layer (mol m^) 

reference oxygen concentration (mol m^?) 
diameter of fiber (m) 

side length of squared duct (m) 

diffusion coefficient of component m (n? s7!) 
oxygen diffusion coefficient (m? s7!) 

effective diffusion coefficient of component m in 
porous media (m? s-1) 

electromotive force (V) 

coefficient in local equilibrium distribution func- 
tion 

velocity distribution function 

local equilibrium distribution function 
Faraday's constant (C mol`!) 

current density (A m?) 

exchange current density (A m^?) 

permeability (m?) 

length and GDL thickness (m) 

step number of random walk 

trial number of random walk 

pressure in LBM 

differential pressure (Pa) 

distance of walker from first position (m) 
mean-square displacement (m?) 

gas constant (m? Pa mol`! K7!) 

resistance of proton transfer through the elec- 
trolyte membrane (2 m2) 

slope of Ap/U calculated by LBM with GDL 
slope of Ap/U calculated by LBM without GDL 
time 

time step 

temperature (K) 

local velocity in LBM 

superficial flow velocity (m s7!) 

voltage 

position 

mesh width 

parameter of Eq. (11) 

transfer coefficient 

variable defined in Eq. (15) (A m mol!) 
porosity 

parameter of Eq. (11) 

viscosity (Pas) 

density in LBM 

tortuosity 

relaxation time 

kinematic viscosity in LBM 


water management become important in regard to its mass trans- 
fer. The proton conductivity of the membrane depends on the 
moisture content, and the proton conductivity rises as the mem- 
brane is wet. Accordingly, supplied gas is usually humidified 
by humidifier, and the humidity is kept high. It is necessary to 
maintain the relative humidity in a cell about 100% in order 
to obtain high output power. However, water condenses when 
relative humidity exceeds 100% and water droplets appear in 
GDL and gas channels. The influence of liquid water on cell 
performance is shown in Fig. 1. When there is liquid water, it 
prevents reactant gas from diffusing and flowing to electrode, 
and the mass transfer rate becomes low. As a result, cell volt- 
age drops because of concentration overvoltage. Therefore, it is 
very important to decrease the bad influences by liquid water 
and to know the two-phase condition in GDL. Generally, water- 
repellent GDL is used in order to improve the performance of 
discharging water. The contact angle of droplet on the surface of 
GDL is almost more than 90°, and it has hydrophobicity. How- 
ever, the detailed correlation between the amount of remaining 
liquid water in GDL and the output performance has not been 
grasped yet. So the direct observations or measurements of liquid 
water remaining in GDL have been examined experimentally. In 
order to measure the two-phase flow condition in GDL which 
is microscopic and opaque in “in situ" condition, the advanced 
techniques are required, for example neutron imaging [1—3]. On 
the other hand, the porous structure and the property of flow and 
mass transfer in single GDL have been experimentally examined 
in some studies as a basic evaluation of GDL [4—6]. To put it 
concretely, the pore distribution, the gas or water permeability 
and the tortuosity have been evaluated by experiments of flow 
and mass transfer. As these parameters depend on the kind of 
GDL structure, the diameter of a carbon fiber, the porosity, the 
clamp pressure, the PTFE content, these effects have also been 
examined. 

Recently, the internal phenomena are examined by the numer- 
ical analysis with the flow and mass transfer property obtained 
by experiments. The effective diffusion coefficient is calculated 
by correcting the diffusion coefficient in a free space with poros- 
ity, tortuosity and interstitial water saturation. Accordingly, it is 
needed to measure how much water content is in GDL to know 
the mass transfer rate of reactant gas (hydrogen or oxygen) from 
gas channels to catalyst layers. It is so difficult to know this value 
that some numerical analysis models with two-phase flow con- 
dition in GDL have been proposed [7—10]. The interstitial water 
saturation has been estimated with various operating conditions 
and current density which affect the rate of water generation 
and condensation, and the influence of liquid water on the volt- 
age drop have been examined. In our past study, the influence 
of GDL thickness on current density distribution was examined 
by numerical analysis [11]. And in [12,13], the mass transfer 
and the flow in the gas diffusion layers were calculated, and 
the approximate model for the GDL mass transfer based on the 
theoretical model was developed. Next, the PEFC reaction and 
thermal flow analysis model which enabled us to calculate an 
actual-sized cell was made with this model. The numerical anal- 
ysis made it possible to examine the separator depth, the GDL 
effective porosity, the GDL permeability, and the flow rate of 
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Fig. 1. Influence of liquid water on internal phenomena. 


the cathode gas which effected on the output performance and 
the current density distribution. However, these analysis models 
simplified the void structure of GDL, and used uniform poros- 
ity and permeability. There are few studies which focus on a 
heterogeneous void structure. We thought that the water perme- 
ation in hydrophobic porous media was strongly affected by this 
heterogeneous structure, and that it was important to examine 
this heterogeneous structure. In this study, in order to calculate 
the water permeation in GDL which consists of heterogeneous 
void structure, as the first stage of the study, simulated GDL 
was made by numerical analysis considering the fiber place- 
ment. And pore size distribution, permeability and tortuosity of 
this simulated GDL were calculated and compared with other 
studies. The validity of this structure was examined from the 
viewpoint of the characteristic of flow and mass transfer on a 
scale of carbon fiber. In addition, it is thought that these numer- 
ical evaluations of the structural property are very important to 
examine the internal phenomena in GDL and to design GDL 
structure. Generally, GDL structural properties have been deter- 
mined by basic experiments, and the data of other references 
have been often used. In this case, it is not easy to understand 
the effect of the size of carbon fiber and the lamination in detail. 
So the effect of the diameter of carbon fiber and the porosity on 
GDL structural properties was examined by numerical analysis. 
In addition, effective diffusion coefficient was calculated from 
tortuosity in simulated GDL, and PEFC output performance was 
evaluated by a simple model. The results of two-phase flow 
numerical analysis based on this simulated GDL will be reported 
in our next paper. 


2. Development of simulated gas diffusion layer 
2.1. Gas diffusion layer in this study 


GDL should be porous media and have high electrical 
conductivity. In general, GDL consists of fiber materials; car- 
bon paper (non-woven fabric) and carbon cloth. These carbon 
materials are distinguished by the difference of the laminated 
structure of the fibers. GDL has water repellent by impregnation 
of PTFE in order to improve the drainage function. In this study, 
carbon paper was examined. It is composed of straight carbon 
fibers dispersed like a plane [14], so the void structure is very 
complex and not uniform. 


2.2. Hypotheses 
The simulated GDL was made on the assumptions as follows: 


1. Carbon fibers are straight and laid in a plane, and they are 
not oriented to the direction of GDL thickness. 

2. Length of carbon fibers is infinite in the target domain, and 
fibers do not have cut surfaces. 

3. Carbon fibers penetrate mutually when they intersect with 
each other. 

4. The diameter of carbon fibers is uniform and constant in the 
same GDL. 

5. Compression and deformation of fibers by clamp pressure 

are ignored. 

. Microporous layer (MPL) is ignored. 

7. Deformation of a porous structure by water repellent with 
PTFE is ignored. 


On 


Actual carbon fibers in real GDL are not in horizontal planes, 
and some fibers are oriented to the direction of GDL thickness. 
However, in this study, it was assumed that the influence of fibers 
oriented to the direction of GDL thickness on the characteristic 
of flow and mass transfer was very low and that most of the fibers 
were set in a two-dimensional plane because of clamp pressure 
when single cell was assembled. And the development of simu- 
lated GDL was simplified. Moreover, the influence of MPL and 
PTFE content, which are important parts of GDL property, were 
ignored. These influences will be examined in our next study. 


2.3. Development of simulated GDL 


On the basis of above assumptions, the simulated GDL was 
developed by a numerical analysis. The setting parameters were 
thickness, porosity, diameter of a fiber and degree of orientation, 
and those values are shown in Table 1. In this study, square area 
with sides 200 um was targeted on the horizontal plane, because 
it was enough to confirm the reproducibility of structural prop- 
erties, as pore size distribution, which would be explained in the 
next section. In addition, the simulated GDL with square 400 jum 
on a side was also made in the base condition shown in Table 1. 
And the simulated GDL was made three times with each setting 
parameters, and the reproducibility was confirmed. First, fibers 
were set horizontally at arbitrary positions and angles in the tar- 
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Table 1 

Setting parameter of development of simulated GDL 

Diameter of fiber (jum) 5,7,9 
Porosity 0.7, 0.6, 0.5 
Thickness (wm) 200 

Size (um?) 200 x 200 
Degree of orientation (°) 360 


The underlined values are basic condition. 


geted domain. The fibers were not cut off in this domain. The 
position and the angle of fibers were obtained from uniform ran- 
dom numbers. This process was continued until the porosity of 
one layer became a setting value in each layer. These fiber layers 
were piled until the GDL thickness became a setting value. By 
the above-mentioned method, the simulated GDL was made. In 
this study, the void part and the solid part were set in a cubic lat- 
tice 1 jum on a side from the viewpoint of the limits of memory 
capacity, computational speed and imaging speed. Accordingly, 
itis necessary to notice that fibers have a ragged surface micro- 
scopically. Fig. 2 shows the simulated GDL in a stereoscopic 


x [um] 


0 100 200 300 400 
Fig. 2. Structure of simulated GDL developed by numerical analysis in the base 
condition: (a) stereoscopic figure, (b) plane figure. 


figure and a plane figure in the case of base condition. In this 
figure, this method numerically developed the porous structure 
which was similar to real GDL. Fig. 3 shows these structures with 
each porosity and fiber diameter. In this figure, it was found that 
the size and the heterogeneity of pore distribution were different 
from those in other conditions. As stated above, although it was 
confirmed qualitatively and visually that the simulated GDL was 
almost same as a real structure observed by a microscope, it was 
needed to evaluate the porous property and to be compared with 
that in other studies in order to confirm the validity of the struc- 
ture. In the next section, the evaluation of the structure was tried. 


3. Evaluation of simulated gas diffusion layer 


As it is difficult to compare directly the microscopic void 
structure of simulated GDL with that of real GDL, macro- 
scopic properties were evaluated. From the point of view of 
two-phase flow and mass transfer, porosity distribution, surface 
roughness, pore size distribution, permeability and tortuosity of 
simulated GDL, which was developed in the previous section, 
were obtained by the numerical analysis. These results were 
compared with those of other studies. 


3.1. Porosity distribution 


In-plane local porosity distribution was examined in order to 
confirm the reproducibility of simulated GDL. This local poros- 
ity was calculated by following method: the average porosity 
between both GDL surfaces in the direction of thickness was 
calculated at 1 jum square mesh which was the lattice size in this 
calculation, and it was obtained at all mesh on GDL surface. 
Fig. 4 shows the local porosity distribution in the base condi- 
tion, and Fig. 5 shows the relationship between the statistics of 
in-plane local porosity and the average porosity with fibers in 
various diameters. It was found that local porosity was almost 
homogeneously distributed with a certain degree of error from 
Fig. 4 and that its statistics became the Gaussian distribution 
from Fig. 5. So, fibers did not unevenly distribute, and they 
were dispersed enough. And it was confirmed that the simulated 
GDL developed in this study was large enough to reproduce the 
macroscopically averaged structure, and that it was possible to 
calculate and compare the various structural properties. On the 
other hand, in Fig. 5, it was found that unevenness of local poros- 
ity was reduced by decreasing the fiber diameter. Thus degree 
of inhomogeneity of average porosity could be evaluated. The 
other cases showed similar results. 


3.2. Surface roughness 


Next, surface roughness of simulated GDL was evaluated. 
It is thought that the surface roughness closely relates to the 
performance of discharging liquid water in gas channels by 
hydrodynamical drag of gas flow. In this study, this property 
was calculated by the local depth from a GDL surface. Fig. 6 
shows the local depth distribution in the base condition. In this 
figure, it was found that some positions are locally deeper owing 
to the condition of lamination. Fig. 7 shows the relationship 
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Diameter of fiber 


Porosity 
5 um 


7 um 


100 — 150. 
x[um] 


Fig. 3. Structures with each porosity and fiber diameter. 


Local porosity [-] 


Fig. 4. Local porosity distribution in the base condition; fiber diameter 7 um, 
average porosity 0.7, thickness 200 jum. 


between the statistics of local depth and the average porosity 
with fibers in various diameters. In this graph, local depth of 
a surface was changed by the diameter of a fiber and the aver- 
age porosity, and the ratio of depth could be expressed by one 
curve. In addition, it was confirmed that maximum depth was 
almost equal to thickness of GDL in the case of porosity 0.7, 
fiber diameter 9 wm, and that some parts between both sides 
of GDL could directly penetrate. Meanwhile, it is needed to 
consider the depth of this surface from the other points. The 
deeper parts become the paths which are straight toward another 
side for flow and diffusion, and it is thought that such parts 
are apparently effective in the mass transfer. However, there 
is liquid water inside GDL when power is generated. Because 
capillary pressure is very large (capillary number is very small), 
the transport of liquid water and mass is more strongly affected 
by local pore size than by the straightness of local path on the 
basis of facts that capillary pressure is inversely related to pore 
radius. This pore size distribution is described in the next sec- 
tion, and the two-phase condition in GDL will be examined 
in the next paper. As there are few other studies about surface 
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Fig. 5. Relationship between the statistics of in-plane local porosity and the average porosity with three sizes of fiber diameters: (a) 5 jum, (b) 7 wm, (c) 9 wm. 


roughness, comparison and verification were not carried out in 
this study. In our future study, it is needed to confirm the validity 
by experiments. 


3.3. Pore size distribution 


The pore size distribution of GDL has been recently dis- 
cussed in other studies about the structural evaluation and the 
behavior of two-phase flow, and it is measured experimentally, 
for example mercury porosimetry. This value is an important 
factor for structural evaluation and is used to examine the per- 
meation of gas or liquid water. In this study, it was obtained by 
calculations as follows: the domain which was enclosed with 
carbon fibers on the plane was extracted, and the area of this 
domain was obtained; next, the equivalent diameter of the area 
was calculated, and the calculated value was regarded as a pore 
diameter. These methods were shown in Fig. 8. In this way, 


the pore size distribution was obtained by statistics of com- 
position ratio of these pore diameters of all closed domains. 
Fig. 9 shows the pore distribution of a simulated GDL in vari- 
ous structural conditions. In these graphs, it was confirmed that 
the structure of the simulated GDL had reproducibility, because 
the pore distributions of the simulated GDL which were made 
three times under the same condition almost agreed mutually. 
The pore size distributions have one peak, and the values were 
about 10-25 um. These values were almost equal to actual mea- 
surements by the other studies [15,16]. Concerning this result, it 
has to be noticed that the two-dimensional pore size distribution 
obtained by void areas in a plane can not be compared sim- 
ply with the three-dimensional pore size distribution obtained 
by experiments, because the pore diameter measured by mer- 
cury porosimetry with Washburn's equation is affected by the 
influence of a steric structure. However, it was found that its 
distribution obtained by this study's method was close to the 
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Fig. 6. Local depth distribution in the base condition; fiber diameter 7 jum, 
average porosity 0.7, thickness 200 jum. 


experimental data in the case of anisotropic porous media such 
as non-woven fabric. In addition, it was found that unevenness 
of pore was reduced by decreasing the fiber diameter and the 
average porosity. 


3.4. Permeability 


Permeability is a parameter of gas flow in porous media, and 
it has been evaluated as the characteristic of GDL by flow experi- 
ments. In this study, the flow analysis in this simulated GDL was 
carried out by lattice Boltzmann method (LBM), and the perme- 
ability was calculated by the relationship between flow rate and 
differential pressure. LBM is a numerical fluid analysis based 
on statistical mechanics [17,18]. In this method, it is assumed 
that continuum fluid is aggregate of numbers of virtual particles 
which moves only to some definite directions. The collision and 
the translation of each particle are calculated with the velocity 
distribution function, and the macroscopic flow field (pressure 
and velocity) is calculated with the moment of this function. 
Therefore LBM has analogy with a kinetic theory of gas which 
treats fluid in the microscopic position. The characteristics of 
LBM are shown as follows: 


e As the collisions and the translations are only repeated, the 
algorithm is simpler than other methods. 

e Unlike a finite differential method and a finite element 
method, equally spaced mesh is used in LBM, so the grid 
generation is not needed. In addition, setting a boundary con- 
dition of walls simpler, LBM is suitable for calculations in a 
complex structure. 

e Unlike a finite differential method and a finite element 
method, the calculation is faster because pressure can be 
explicitly calculated by density. 


Moreover, the complex shape of an interface can be cal- 
culated, so this method has been applied to a two-phase flow 
[19-21]. 

In this study, three-dimensional 15 velocities model (3D15V 
model), which is shown in Fig. 10, was used as lattice models 
of LBM. And BGK model was used as the model of collision 
term. At position x and time f, velocity distribution function 
fm Gt), expresses existing probability of virtual particle moving 
at a velocity Cm, satisfies the following lattice Boltzmann 
equation: 


fm (x + Cm Ax, t + At) — fin(x, t) 


— Sin, t) — fa (x, t) (1) 


Tr 


where fmm! is the local equilibrium distribution function, tr the 
relaxation time, Ax and Aż are the mesh width and the time 
step, respectively, and these values are one in LBM. In this 
equation, the left-hand side means the translation of virtual 
particle, and the right-hand side means the change of velocity 
distribution function by collisions. The local equilibrium 
distribution function is calculated by the following equation: 


(x, t) = Emp |1 + 3€m -u + (en .uy — ju -u (2) 
where u is the local velocity, p the density. The coefficient 
Em is fixed in order to satisfy Navier-Stokes equation and the 
definitional equation of density and velocity, and these values 
of 3D15V model are shown in Table 2. Fluid density o and flow 
velocity u, which are macroscopic variable, are defined with 
velocity distribution function fm as follows: 


N 
P=} fn (3) 
m=1 
1 N 
u- -Y c, f, (4) 
p m=1 
Table 2 
Coefficient E,, of 3D15V model 
m=1 En =2/9 
m=2-7 En = 1/9 
m=8-15 Ey, =1/72 
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Fig. 7. Relationship between the statistics of local depth and the average porosity with three sizes of fiber diameters: (a) 5 jum, (b) 7 wm, (c) 9 jum. 
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Fig. 8. Schematic diagram of method to obtain pore size distribution. 


G. Inoue et al. / Journal of Power Sources 175 (2008) 145—158 153 


(à o 


IS 
E 6 
2 
£ 
g 
o 
a 4 
A 1 
2 ' v SIX 
ull || sr ^d 
© 
E 
lali: 
N 
0 E Y P will Pea 
0.1 10 100 1000 
pore diameter [um] 
| | il 
9 
rj 
£ 
5 
A 


pore diameter [um] 


S 
2 
£ 
S 
2 
S 
KA 
Vy d 
Sul 
* a | 
0.1 1 10 100 1000 
pore diameter [um] 
(d) 
9 
2 
g 
$ 
S 
X n jl 


0.1 1 10 100 1000 
pore diameter [um] 


Fig. 9. Pore size distribution of simulated GDL with various sized diameter and porosity of fibers: fiber diameter (a) and (b) 7 jum, (c) and (d) 9 jum; porosity (a) and 


(c) 0.7, (b) and (d) 0.6. 
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Fig. 10. Lattice structure of three-dimensional 15 velocities model (3D15V 
model) of LBM. 


Kinematic viscosity v is calculated as the function of relaxation 
time. 


E 1 
v=; (7-3) © 


As mentioned above, LBM can be applied to direct fluid 
analysis in porous media with complex structure such as GDL. 
However, as all variables are dimensionless values in LBM, it is 
not possible to calculate them with concrete physical quantity. 
Accordingly, the actual value of pressure drop and permeability 
cannot be calculated directly. So the permeability of simulated 
GDL was obtained by the following method in this study. In 
Fig. 11, simulated GDL was set in a squared duct, and supplied 
gas flowed from the downside in various superficial flow 
velocity (U). Then, differential pressure between both sides 
of GDL (Ap) was obtained by LBM. In calculation by LBM, 
these differential pressure and flow rate are dimensionless 
values. Assuming the gas flow in simulated GDL could be 
approximated by Darcy’s law, the slope SgpL, which equalled to 
Ap/U, was calculated by a least-squares method. Similarly, the 
slope Sref was calculated in the case of no GDL, as reference. 
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Simulated GDL 


"A 


Uniform inlet flow 


Fig. 11. Calculation model of gas flow in GDL for evaluation of permeability. 


On the other hand, the relationship between pressure drop and 
gas flow velocity in porous media is theoretically calculated by 
the following Darcy's law [22]: 


uL 
Ap = —U 7 
P= (7) 


where K is the permeability, u the viscosity and L the section 
length. In addition, the pressure drop in a squared duct without 
porous media is calculated by the following equation which is 
derived from Hagen-Poiseuille law and a friction coefficient of 
the squared duct [23]: 


28.44 uL 
Dis the side length of squared duct. Assuming the ratio between 


SGDL and Sief obtained by LBM was the same that between 

Eqs. (7) and (8), the following equation was derived: 
. uL 28.44uL 

SGDL : Stef = K : C39 


By developing the equation above, the following permeability 
equation is obtained: 


(8) 


(9) 


E D? Sref 
~ 28.44 SGDL 


(10) 


In this study, the in-plane and through-plane permeability of sim- 
ulated GDL was calculated by the process above, and the rela- 
tionship between permeability and porosity was calculated with 
various GDL. Furthermore, this result was compared with that 
of other studies in order to verify the validity. Gostick et al. [24] 


measured the in-plane and through-plane permeability of various 
GDL, and the effect of porosity was examined by compressing. 
Besides, it was confirmed that these results exactly agreed to the 
results of the following Tomadakis-Sotirchos model (TS model) 
which estimated permeability of fibrous porous media [25]: 


B € (e — e) 2a? 
— 32(n e)? (1— epY*[(o + De — £p]? 


(41) 


where dr is the diameter of fiber, € the porosity, œ and £p are the 
constants that depend on the structure and on the flow direction. 
When randomly set fibers form a plane and such planes are 
parallel to each other, œ and £y are respectively 0.521 and 0.11 in 
the case of in-plane flow. And in the case of through-plane flow, 
these values are 0.785 and 0.11, respectively. Fig. 12 shows an 
example of flow analysis results by LBM in a simulated GDL. 
In this graph, flow velocity distribution in porous media which 
is a complicated shape was calculated by this method. These 
analyses were carried out on other GDL, and the differential 
pressure and permeability were obtained. Fig. 13 shows the 
calculation results and TS model standing for a relationship 
between in-plane and through-plane permeability and porosity 
in the case of fiber diameter 5 jum and 9 jum. In these graphs, it 
was found that the calculation results were almost equivalent to 
TS model in the case of in-plane flow, and that validity could be 
confirmed. So, this calculation method was effective to predict 
permeability. However, the calculation results of through-plane 
permeability were different from results obtained from a model 
equation; the former was twice as much as the latter. And the 
permeability in both directions was almost equal to each other. 
Though the permeability of an actual non-woven fabric has 
anisotropy, it was thought that this difference arose owing to the 
insufficiency of calculating area and low-resolution mesh in this 
study. Nonetheless the approximate permeability equivalent to 
model equations could be predicted by this analysis method 
though it had a certain amount of error. 


3.5. Tortuosity 


Tortuosity is an important parameter of mass transfer in 
porous media. The effective diffusion coefficient in porous 
media with porosity and tortuosity is shown as follows: 


€ 
DË = Dm : (12) 


Generally, tortuosity of isotropic porous media, such as spherical 
particle, is obtained by Bruggeman equation [26]. However, this 
equation cannot be used in the case of non-woven fabric because 
of its anisotropic structure. Tomadakis and Sotirchos calculated 
tortuosity of fibrous porous media and modeled its equation by 
random walk method [27]. In this study, the validity of tortuosity 
in simulated GDL was confirmed by random walk method. In the 
calculation by Tomadakis and Sotirchos, tortuosity was obtained 
when the overall direction of mass transfer was given previously. 
In this study, as a supplementary examination, tortuosity was cal- 
culated in the case that the mass transfer direction was not given 
by the method of Watanabe and Nakashima [28]. In simulated 
GDL, walker takes a step of lattice size to the arbitrary direction 
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Fig. 12. Example of flow analysis results by LBM: (a) GDL structure, (b) pressure distribution, (c) flow rate distribution, and (d) velocity distribution. 


determined by random numbers from an arbitrary vacant point. 
Walker does not move to solid phase. Walker moves m steps, and 
the distance r between the first position and the walker’s posi- 
tion at step m is obtained. This process was tried n times, and 
the mean-square displacement (7?) was calculated as follows: 


2 c : 2 
acr (13) 


Similarly, mean-square displacement in free space was cal- 
culated, and the tortuosity was obtained by dividing the value 
in free space by that in simulated GDL with various average 
porosities and diameters of fiber. As the assumption of a method 
of Watanabe and Nakashima [28] was that porous structure had 
isotropy, it should be noticed that the results of this calculation 
were not always equal to the actual tortuosity of GDL with 
anisotropy. Fig. 14 shows the relationship between porosity and 
tortuosity of simulated GDL, and it also shows the Tomadakis 
and Sotirchos model which is the model of fibrous porous 


media [27]. In this graph, it was found that calculation results of 
this study were almost equal to the value between in-plane and 
through-plane results of TS model. As the direction of diffusion 
was not determined in this calculation, these calculated values 
became intermediate value between both values of TS model. 
As a result, the validity of winding structure of this simulated 
GDL was confirmed qualitatively. In addition, there are some 
results which need to be focused on. Tortuosity ought to be 
independent of the diameter of fiber as long as geometrical 
similarity is kept. However, tortuosity was changed by the 
diameter of a fiber in this calculation. This result was caused by 
insufficiency of resolution, and it needs to be calculated with 
high resolution in the next study. 


3.6. Evaluation of cell output performance with calculated 
tortuosity 


PEFC cell output performance was evaluated by simple 
model including calculated tortuosity in simulated GDL. The 
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Fig. 13. Calculation results and TS model of relationship between in-plane and 
through-plane permeability and porosity with two sizes of fiber diameter: 5 jum 
(a), and 9 jum (b). 
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Fig. 14. Calculation results and TS model of relationship between tortuosity 
and porosity with various sized fiber diameter. 


model equations were derived from our past model [12]. The 
oxygen diffusivity in cathode GDL was modeled, and the volt- 
age was calculated by various overvoltage. In this study, it was 
assumed that there was only gas phase and that temperature was 
constant and uniform. In addition, the electrolyte membrane was 
humidified well, and ionic conductivity was constant and uni- 
form. Though our past model included the gas flow through 
GDL, this flow was ignored in this study. The overvoltage of 
anode reaction was ignored, and the relationship between cell 
voltage and overvoltage is shown by the following equation: 


RT i 
n 
2 F io, 0» 


< 
Il 
t3 


Rohmi (14) 


where F is the electromotive force, R the gas constant, T the tem- 
perature, a; the transfer coefficient, F the Faraday’s constant, i 
the current density, io, the exchange current density, C or the ref- 
erence oxygen concentration, Co, the oxygen concentration at 
an interface of catalyst layer and Rohm is the resistance of proton 
transfer through the electrolyte membrane. Rohm was calculated 
by Nguyen’s equation [29]. Reference oxygen concentration and 
exchange current density were shown as the variable y by the 
following equation: 


_ io 
Y= Ca (15) 
The oxygen diffusion in GDL was simplified with the following 
equation which was derived from linear concentration gradient 
in Fick's law: 


Cò, = M coanel — MEO ERR S (16) 


where coe is the oxygen concentration at a channel, L the 
thickness of GDL, Do, the oxygen diffusion coefficient, t and € 
are tortuosity and porosity in GDL. These tortuosity and poros- 
ity were obtained by the evaluation of simulated GDL in the 
previous section, and the effect of GDL structure on output 
performance was examined. 

Calculation condition was shown in Table 3. However, as this 
model did not include the effect of liquid water, multi-content 
diffusion and heat transfer, it was not possible to compare the 
absolute values. So the maximum current density when porosity 
was 0.7 with simulated GDL was used as the reference value, and 
the output performance of each GDL was compared by using rel- 
ative current density. Fig. 15 shows the current density-voltage 
curve of each GDL arranged by relative current density in the 


Table 3 

Calculation condition of output performance 

Temperature (°C) 60 

Transfer coefficient 0.3 

Porosity 0.5-0.7 
Electromotive force (V) 1.23 

GDL thickness (um) 200 
Membrane thickness (jum) 30 

Oxygen diffusion coefficient (m? s~!) 2.00 x 10? 
y of Eq. (15) (Am mol!) 4.0 x 10? 
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Fig. 15. Relative current density—voltage curve calculated from each tortuosity 
in the case of various porosities. 


case of various porosities. Fine lines show the calculation results 
with tortuosity in simulated GDL in the case of 5 jum-diameter 
fiber, and bold lines show the Bruggeman equation. In this figure, 
it was found that the tortuosity affected to the output property, in 
particular, the influence was stronger in the case of low poros- 
ity. It was thought that the calculated results with Bruggeman 
equation which was generally used in other papers overestimate 
the diffusivity of GDL which consists of fibrous porous media. 
Moreover, though this model did not include the effect of liquid 
water, it is expected that it is important to know the substantial 
tortuosity relating the diffusivity in the case of two-phase condi- 
tion. In this study, this model was so simple that the quantitative 
comparison could not be carried out. It will be needed to develop 
the high-accuracy model. 


4. Conclusion 


As a basic study for examination of flow and mass trans- 
fer in GDL with heterogeneous porous structure, the simulated 
GDL was developed by numerical analysis. The parameters 
were the diameter of a fiber, the average porosity and thickness. 
Furthermore, the structural properties, such as local porosity 
distribution, surface roughness, pore size distribution, perme- 
ability, tortuosity, were evaluated, and the validity was confirmed 
by comparing with the results of other studies. As a result, the 
effects of average porosity and the diameter of a fiber on these 
structural properties were clarified by numerical analysis. The 
prediction method of this study can be applied to other porous 
structure except for non-woven fabric, and it is thought that 
new GDL structure which contributes to the high power density 
and the high durability can be developed and proposed by this 
model. Moreover, optimal GDL structure can be examined by 
applying this method from the viewpoint of gas diffusivity and 
water discharging efficiency in GDL. Furthermore, the effective 
diffusion coefficient was calculated from tortuosity in simulated 
GDL, and PEFC output performance was evaluated with a sim- 
ple model. It was confirmed that the diffusion overvoltage was 
largely affected by tortuosity on a general operating condition 
and that it was important to calculate the tortuosity with accu- 


racy. However, as this evaluation of output performance did not 
include the effect of liquid water, multi-content diffusion and 
heat transfer, it will be needed to develop the model including 
these phenomena with simulated GDL in our future study. In 
addition, as GDL was transformed by clamp pressure, the influ- 
ence of transformation will be considered. On the other hand, 
though pore diameter distribution was obtained as planate pore in 
this study, it is expected that pore size has three-dimensionality 
in à cubic structure. In our future study, it needs to be improved 
the prediction method of pore diameter distribution. 


Acknowledgement 


This research was supported by the research and develop- 
ment of polymer electrolyte fuel cell from the New Energy 
and Industrial Technology Development Organization (NEDO), 
Japan. 


References 


[1] D. Kramer, J. Zhang, R. Shimoi, E. Lehmann, A. Wokaun, K. Shinohara, 
G.G. Scherer, Electrochim. Acta 50 (2005) 2603-2614. 

[2] J. Zhang, D. Kramer, R. Shimoi, Y. Ono, E. Lehmann, A. Wokaun, K. 
Shinohara, G.G. Scherer, Electrochim. Acta 51 (2006) 2715-2727. 

[3] M. Hickner, N. Siegel, K. Chen, D. Hussey, D. Jacobson, M. Arif, Fuel Cell 
Seminar, Honolulu HI, November 13-17, 2006, pp. 380—383 (Abstracts). 

[4] J. Ihonen, M. Mikkola, G. Lindbergh, J. Electrochem. Soc. 151 (2004) 
A1152-A1161. 

[5] M.V. Williams, E. Begg, L. Bonville, H.R. Kunz, J.M. Fenton, J. Elec- 
trochem. Soc. 151 (2004) A1173—A1180. 

[6] J.G. Pharoah, J. Power Sources 144 (2005) 77-82. 

[7] U. Pasaogullari, C.Y. Wang, K.S. Chen, J. Electrochem. Soc. 152 (2005) 
A1574-A1582. 

[8] U. Pasaogullari, C.Y. Wang, J. Electrochem. Soc. 151 (2004) A399-A406. 

[9] G. Lin, W. He, T.V. Nguyen, J. Electrochem. Soc. 151 (2004) 
A1999-A2006. 

[10] N.P. Siegel, M.W. Ellis, D.J. Nelson, M.R. von Spakovsky, J. Power Sources 
128 (2004) 173-184. 

[11] G. Inoue, Y. Matsukuma, M. Minemoto, J. Power Sources 154 (2006) 
8-17. 

[12] G. Inoue, Y. Matsukuma, M. Minemoto, J. Power Sources 157 (2006) 
136-152. 

[13] G. Inoue, Y. Matsukuma, M. Minemoto, J. Power Sources 157 (2006) 
153-165. 

[14] M.F. Mathias, J. Roth, J. Fleming, W. Lehnert, Handbook of Fuel Cells, 
vol. 3, Wiley, 2003, pp. 517—537 (Chapter 42). 

[15] J.T. Gostick, M.W. Fowler, M.A. Ioannidis, M.D. Pritzker, Y.M. 
Volfkovich, A. Sakars, J. Power Sources 156 (2006) 375—387. 

[16] H.K. Lee, J.H. Park, D.Y. Kim, T.H. Lee, J. Power Sources 131 (2004) 
200-206. 

[17] S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, 
Oxford Science Publications, New York, 2001. 

[18] D.A. Wolf-Gladrow, Lattice-Gas Cellular Automata and Lattice Boltzmann 
Models, Springer, New York, 2000. 

[19] T. Yoshimoto, G. Inoue, Y. Matsukuma, M. Minemoto, Proceeding of 
Renewable Energy, Makuhari Japan, October 9-13, 2006, pp. 1364- 
1367. 

[20] Y. Matsukuma, S. Sakashita, G. Inoue, M. Minemoto, Proceedings of the 
15th Discrete Simulation of Fluid Dynamics (DSFD 2006) Conference, 
2006. 

[21] T. Koido, T. Furusawa, K. Moriyama, K. Takato, 210th Meeting of the 
Electrochemical Society, Cancun, Mexico, October 29-November 3, 2006. 

[22] D.A. Nield, A. Bejan, Convection in Porous Media, Springer, New York, 
1998. 


158 G. Inoue et al. / Journal of Power Sources 175 (2008) 145—158 


[23] R.B. Bird, W. Steward, E.N. Lightfood, Transport Phenomena, Wiley, New [26] J. van Brakel, P.M. Heertjes, Int. J. Heat Mass Transfer 17 (1974) 


York, 1960. 1093-1103. 
[24] J.T. Gostick, M.W. Fowler, M.D. Pritzker, M.A. Ioannidis, L.M. Behra, J. [27] M.M. Tomadakis, S.V. Sotirchos, AIChE J. 39 (1993) 397-412. 
Power Sources 162 (2006) 228-238. [28] Y. Watanabe, Y. Nakashima, Comput. Geosci. 28 (2002) 583-586. 


[25] M.M. Tomadakis, T.J. Robertson, J. Compos. Mater. 39 (2005) 163-188. [29] T.V. Nguyen, R.E. White, J. Electrochem. Soc. 140 (1993) 2178-2186. 


